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Abstract 

The development of multiple-relaxation-time (MRT) Lattice Boltzmann method (LBM) is a sig¬ 
nificant contribution in improving the numerical behavior, revealing the math and physics mecha¬ 
nism and extending the application of LBM. However, some of the MRT schemes proposed previ¬ 
ously are not physically-consistent. In this work, we take D2Q9 as a example to show how to derive 
physically-consistent MRT-LBM schemes by eigenvalue decomposition of the collision operator. In 
addition, the scheme is validated by the equivalence to Navier-Stokes equations and numerical 
simulations. 
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I. INTRODUCTION 


The past two decades have seen the rapid growth of Lattice Boltzmann method (LBM) |T] pfj. 
Among the many contributions, the proposal of multiple-relaxation-time (MRT) collision 
model takes its place. The significance of MRT-LBM is threefold: it improved the numerical 
behavior of LBM [3], reveals the math and physics behind LBM a 0. and facilitate the 
extension of LBM j6] [7] [8]. 

Roughly speaking, the MRT collision model is an extension of the BGK collision model 
by decomposing the collision process into different modes and assigning different parameters 
for each mode. Obviously, the modes with different relaxation times should be independent, 
otherwise, inconsistency will happen. However, this rule is not always obeyed in practice. In 
this work, we will take D2Q9 lattice as an example to show how to get physically-consistent 
MRT schemes. 

The rest of the paper will be organized as follows: in Section [TTJ a brief introduction on 
LBM will be given with an emphasis on MRT collision model; in Section IH we will propose 
a physically-consistent MRT scheme based on an eigenvalue decomposition on the linear 


approximation of the BGK collision operator; in Section |IV[ we will prove that the scheme 
reduces to Navier-Stokes equations at the macroscopic level; and the scheme will be further 
validated by simulation in Section [Vj finally, we will conclude the work by Section VI 


II. LATTICE BOLTZMANN METHOD 

In this section, a brief introduction on Lattice Boltzmann method (LBM) will be given. 
Throughout the paper, scalars, vectors and tensors are denoted by lowercase letters, lower¬ 
case letters in boldface and uppercase letters, respectively. 

In LBM, the flow of fluids is simulated by particles hopping on a lattice. Fixing the time 
step At and the unit length of the lattice Ax, the velocity of the particles can only be chosen 
from a finite set of vectors {e* | i — 1,..., N}. Let f(x, t ) be the particle velocity distribution 
function whose i-th component /)(x, t) gives the portion of particles with velocity e* at node 
x at time t. It satisfies the Lattice Boltzmann Equation (LBE) 

/i(x + ejAi,t +At) = /i(x,t) + n, ; (f(x, t)) (1) 

where i = 1,2, ...,1V and n,j(f(x,t)) is the 1-th component of the collision operator that 
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FIG. 1. D2Q9 Scheme 

describe the effect of collision of particles. In coding, the computation of (JT|) is usually 
divided into two steps: 

• Collision: 

/iOM) = /*0M) + n l (f(x,t)), ( 2 ) 

• Streaming: 

/;(x + e;At,t + At) = /-(x,t). (3) 

In this work, we will take the D2Q9 lattice (Figure [T]) as an example, where the set of 
velocity is taken to be 

e 0 = [0, 0], ei = [ k , 0], 

e 2 = [0, k ] , e 3 = [- k , 0], 
e 4 = [0, 0], e 5 = [«, k ] , 
e 6 = [-k, k ] , e 7 = [-«, - k ] , 
e 8 = [k, -k] • 

where ac = Ax/ At is the characteristic velocity of the lattice. The arguments in the rest of 
the paper extend easily to other kinds of lattices. 

The particle velocity distribution function f(x, t) is related to the macroscopic variables 
by 

N N 

p='52fu p u = ( 4 ) 

2—0 2—0 

where p and u are the macroscopic density and velocity respectively. Due to the conservation 
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of mass the momentum, for any f(x, t), the collision operator satisfies 


N 


N 


5Z n «(f(x,t)) = o, J>(f(M))e, = o. 


i =0 


i =0 


( 5 ) 


In general, the particle velocity distribution function f(x, t) cannot be deduced from 
the macroscopic quantities p and u. But in equilibrium, the particle velocity distribution 
function f eq (x, t) only depends on p and u. In the D2Q9 lattice, by fitting f(x, t) to the 
continuous equilibrium distribution function, we obtain that [9] 

/r = P w i ^1 + 3u • ej + ^ (u • e*) 2 - ^|u| 2 ^ (6) 


with wq = 4/9, wi = W 3 = w$ = wj = 1/9 and w\ = W 3 = w$ = W 7 = 1/36. 

Generally speaking, the collision operator is trying to restore the particle velocity distri¬ 
bution function f(x, t) to its equilibrium distribution f eq (x, t). In the BGK collision model, 
the collision operator is taken to be a linear relaxation operator 


n(f(x,t)) 


f eq (x,t) -f(x,f) 
TbGK 


(7) 


where the relaxation time tbgk is determined by kinematic viscosity u by 

1 ( 3z/A t 
T= 2 + ~(KxY' 


( 8 ) 


In the multiple-relaxation-time collision model, the particle velocity distribution function 
f(x, t) is decomposed into N (the dimension of f(x, t)) modes and different relaxation times 
are assigned to each mode, i.e. 


n(f(x,f)) = f>- 1 r- I p(f*>(x,«) - f(x, t)), 


(9) 


where P is an invertible matrix and T = diag (ti, t 2 , ..., Tn). When Ti — ... — Tn — Tbgk , 
the MRT model degenerates to the BGK model. 

The decomposition of f(x, t) is not arbitrary; the modes corresponding to different relax¬ 
ation times should be independent. Previously, the most common decomposition scheme is 
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given by [5] 


P = 


1 

1 

1 

1 

1 

1 

1 

1 

1 

-4 

-1 

-1 

-1 

-1 

2 2 

2 

2 

4 

2 

2 

2 

2 

1 

1 

1 

1 

0 

1 

0 

-1 

0 

1 -1 

-1 1 

0 

-2 

0 

2 

0 

1 -1 

-1 1 

0 

0 

1 

0 

-1 

1 

1 

-1 -1 

0 

0 

-2 

0 

2 

1 

1 

-1 -1 

0 

1 

-1 

1 

-1 

0 0 

0 

0 

0 

0 

0 

0 

0 

1 -1 

1 

-1 


( 10 ) 


where the rows correspond respectively to density, energy, energy square, (E-momentum, x- 
energy flux, ^-momentum, y-energy flux, diagonal component and off-diagonal component of 
stress tensor. However, the modes in this decomposition are not independent, e.g. “energy” 
and “energy square”. Thus, inconsistency may happen when different relaxation times are 
assigned to them. 


III. EIGENVALUE DECOMPOSITION 

Before proposing a physically-consistent MRT scheme, we will first perform an eigenvalue 
analysis on the BGK collision operator. The LBM holds when the macroscopic velocity 
|u(x,f)| <C Ax/At, i.e. the Mach number of the flow with respect to the characteristic 
velocity k is small. Therefore, ([6]) is well approximated by its linearization 

,/f 1 = pWi (1 + 3u • e t ). (11) 


In this case, f eq depends linearly on f by 


f = T 'f eq 


( 12 ) 
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The set of eigenvalues of matrix T are {1,1,1, 0, 0, 0, 0, 0,0} and the corresponding eigen¬ 
vectors are taken to be 


c 1 = [1,1,1,1,1,1,1,1,1], 
c 2 = [0,1,0, -1,0,1,-1, -1,1], 

c 3 = [0,0,1,0, -1,1,1, -1,-1], 

1 2 12 1 2 2 2 2 ' 
C4 _ 3’3’ 3’3’ 3’3’3’3’3 ’ 

c 5 = [0,0,0,0,0,1,-1,1,-1]. 

'1 12 12222 2 ' 
Cg 3’3’ 3’3’ 3’ 3’ 3’ 3’ 3 ’ 

c 7 = [0,1,0,-1, 0,-2, 2, 2,-2], 

c 8 = [0,0,1,0, -1,-2, -2,2, 2], 

c 9 = [-1,0, 0,0, 0,4,4,4,4]. 


(14) 


The nine eigenvectors, each representing a physical mode, are divided into three groups: 
{ci,c 2 ,c 3 }, {c 4 ,c 5 ,c 6 } and {c 7 ,c 8 ,c 9 }. 

In the first group, c 1? c 2 ,c 3 correspond to density, x-momentum and y -momentum re¬ 
spectively, 


p(x,t) = c 3 • f (x,t) 

p(x, t)u x (x, t) = c 2 • f (x, t) (15) 

p(x,t)u y (x,t) = C 3 • f (x, f) 
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By (J2]) and (12), for each k = 1 , 2, 3 and any r k , we have 

c k ■ f'(x, t) = c k • f(x,t). 


( 16 ) 


Therefore, mass and momentum are conservative in the collision step. 

In the second group, c 4 , c 5 , c 6 correspond to the components of the stress tensor T EU by 


Jxx(-x,t)/v = c 4 • f (x,t) 

7xy(x, t)/^ = C 5 • f (x, f) (17) 

7 ra (x,t)/n = c 6 • f(x,t) 


where 


r = 


'T## ^/xy 


(18) 


Ixy 7TO 

The physical picture is explained in the following way. Take the control volume as shown by 
the dashed line in Figure[lJ the flux of r-momentum in the x - dir ection crossing the boundary 
of T is given by 


# ra = [0,l,0,l,0,l,l,l,l]-f(x,i) 


(19) 


Similarly, the flux of y-momentum in the y-direction crossing the boundary of T is given by 


% y = [0,0,1,0,1,1,1,1,1] -f(x,f). (20) 

The flux of a;-momentum in the y-direction crossing the boundary of T is equal to flux of 
y-momentum in the (E-direction, that is, 


= iyx = (0,0, 0,0, 1,-1, 1,-1] ■ f(x,f). 


( 21 ) 


Noting that 


^xx ^xy 

$ $ 

'Vyx ^yy 


= T/u+pi 


( 22 ) 


and p = |, we obtain ([lT]) . 

As shown above, the three eigenvalues in this group correspond to the transportation of 
momentum due to viscosity. Therefore, the relaxation times should be taken as r 4 = 75 = 
re = tbgk ■ By 0 and ([12]), for each k = 4, 5, 6 , we have 


c k ■ f'(x, t) = (1-) c k ■ f(x, t), 


t BGK 


(23) 
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in the collision step. 


Finally, in the third group, c 7 ,c 8 ,c 9 represent no macroscopic physical quantity, namely, 
they are redundant degrees of freedom in computation. Therefore, the relaxation times 
77 , ts, tq can be chosen arbitrarily. By ([2]) and (12), for each k — 7, 8, 9, we have 


c k ■ f'(x,t) = (1-) Cfc • f(x,t), 

1~k 


(24) 


in the collision step. 

Based on the eigenvalue decomposition, we propose the following MRT collision operator 


n(f(x,t)) = -Q~ l SQ f(x,t), 


where 
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The relaxation time tbgk is determined by kinetic viscosity v by 


Obviously, 


1 | 3uAt 

2 + (Ax) 2 


c fc ■ Q X SQ 
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o, 

if k = 1,2,3 

C k/ T BGK, 

if k = 4,5,6 

Ck/l~ki 

if k = 7, 8, 9 
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(25) 


(26) 


(27) 


(28) 
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IV. MULTIPLE-RELAXATION-TIME COLLISION OPERATOR 


The MRT scheme proposed in Section III is validated by the equivalence to Navier-Stokes 
equations through Chapman-Enskog expansion ra- Assuming that Ax — At — e, then the 
Taylor expansion of LBE ([!]) gives 


£ [ + e i ' V fi) + £2 ( o 


2 dt 2 +Gi ' V dt + 2 e « e i - vy /i ) 


(30) 


Let ti be the convection time scale and t 2 be the diffusion time scale, then the time derivative 
decomposes to 

d d d 

dt = dt l +£ dt 2 ' ^ 

Accordingly, the distribution function decomposes near the equilibrium to 


/i(x,i) = /° q (x,i) +ef, m (x,t) ±e 2 ff 2, (x,t) + .... 


(32) 


By conservation of mass and momentum Q, we have 

N N 

= p' ^ ei = pn - 

i=0 i =0 


(33) 


and for s — 1, 2,... 






i=0 


Plugging (31) (32) into (30) and noting that 


(34) 


Q-'.SQ r* = 0, 


(35) 


we obtain that 


to the order e, and 


df, 


dt\ 


eq 

^ + e,.V/r = -[Q- 1 SQf' 1 >] i 


(36) 


a/° g 9fP 

dt 2 dti 
to the order e 2 . Using 


+ e, • V/<‘> + + e,. + re,e, : VV/’ 1 ’ = - [Q-'SQ f <2 >], (37) 


1 a 2 /, 


(1) 


.a// 1 ’ i 


(l) 


2 dt 2 ' ~ l ' dt 
simplifies to 


3^ + dAAAA + 6j . v [ M fq) = - [g-‘so f 2 >], 


dt 2 
-1 


(38) 


where M — I — Q 1 SQ/2. 
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Using (33) (34), the mass equation is given by YliL o ((36) + (38)) as 


! + V.(„u) = 0 


( 39 ) 


and the momentum equation is given by Y^i=o e i' ((36) + (38)) as 

N 


d(pu) 

dt 


+ V-^(e i e i /? q + e i e i [M f (1) ]J=0. 


where 


i =0 


N 


= pl + puu. 


(40) 


(41) 


i=0 


and by (29) 


N 


i M f(i) L = 1 


i =0 


2 tbgk 


C 4 • f ^ C 5 • fW 

C 5 • f (1 ) Cg • f (1) 


= r 


(42) 


From the discussion above, we make the following remarks. Let <S 4 , <S 2 ,<S 3 be the linear 
subspaces spanned by {ci, c 2 , c 3 }, {c 4 , c 5 , c 6 } and {c 4 , c 5 , c 6 , c 7 , c 8 , c 9 } respectively, then 

• ci, c 2 , c 3 is arbitrary as long as <Si is preserved; 

• c 4 , c 5 , c 6 is arbitrary as long as S 2 is preserved; 

• c 7 , c 8 , c 9 is arbitrary as long as S 3 is preserved; 

Though the choice of r 7 , r 8 , r 9 has no influence on the result, for stability, we require that 


U,t 8 ,t 9 > 


(43) 


V. SIMULATIONS 


In this section, the simulation results of 2-D cavity flow will be given to verify the MRT 
scheme proposed in Section IV As shown in Figure [2| the size of the domain is 50 x 50. The 
west boundary x — 0, south boundary y — 0 and east boundary x = 50 are solid walls. The 
north boundary y = 50 moves at a constant speed of (0.1, 0). 

In the simulations, the lattice has 51 x 51 nodes where the unit length Ax = 1 and the 
time step At — 1. Initially, the density and velocity on each node are set to be p = 1 
and (u, v ) = (0, 0) respectively and the velocity boundary conditions are applied at the 
boundary HU- The relaxation times are assigned differently in the following seven cases 
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FIG. 2. Cavity flow 


Case 1 : S — diag (o, 0, 0, —, —, —, —, —, —^ 

^ V T a Ta Ta T a' T 0 ’ T a ) 


• Case 2 : S 


diag ( 0 , 0 , 0 , 
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• Case 6 : S 


diag ( 0 , 0 , 0 , 
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• Case 7: S 


diag ( 0 , 0 , 0 , 
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where u a = 0 . 2 , zy, = 0 . 6 , r a = \ + 3z/j and n = ^ + 3z/ 2 . 

In Figure [3] [ 8 j the plots of ^-velocity against y along the line x = 25 (shown by the 

thick solid line in Figure [2]) at time t = 125 and t — 175 are compared between different 
assignments of relaxation time. The results show that changing the relaxation time of mode 
C 4 , C 5 and Cg affects the computation results, while changing the relaxation time of mode 
C 7 , c 8 and c 9 does not. In addition, the relaxation time of mode c 5 has a greater influence 
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FIG. 3. Change the relaxation time of mode C 4 
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FIG. 4. Change the relaxation time of mode C 5 
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FIG. 5. Change the relaxation time of mode c$ 
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FIG. 6. Change the relaxation time of mode C7 
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FIG. 7. Change the relaxation time of mode cs 



y 


FIG. 8. Change the relaxation time of mode eg 
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than the relaxation time of mode C 4 and Cq. This is because the shear stress plays a more 
important role than the normal stresses in this physical process. 


VI. CONCLUSION 

In this work, we proposed a way of deriving physically-consistent MRT-LBM schemes 
based on eigenvalue decomposition of the collision operator. We showed that the scheme is 
equivalent to the Navier-Stokes equations at the macroscopic level and is in agreement with 
the simulation results. 
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